%  Copyright (C) 2002 Regents of the University of Michigan, portions used with permission 
%  For more information, see http://csem.engin.umich.edu/tools/swmf
\section{Divergence of B Control \label{section:divb}}

There is a big difference between the view of theorists,
who would generally insist that $\divB$ should be exactly zero, and
practitioners of numerical MHD, who usually take a more pragmatic approach,
and are satisfied with $\divB$ converging to zero as
the grid resolution $\Delta x$ and the time step $\Delta t$ approach zero.
The justification for the latter approach is simple:
none of the numerical values agree to the analytical solution exactly,
so why should one insist that a specific combination of them,
namely some numerical representation of $\divB$ should be equal
to the analytic value, i.e. zero.

One way of ensuring a small numerical value for $\divB$ is
to demand that some particular discretization is exactly zero.
The projection and the constrained transport schemes fall into this 
class.
Another possibility is to set the numerical value of $\divB$ to zero in the
initial and boundary conditions, and to trust the scheme to maintain this
condition until the end of the simulation to the accuracy of the truncation
error. The 8-wave scheme and the diffusive divergence $B$ control belong
to this class.

\subsection{8-wave Scheme \label{section:8wave}}

The default scheme used in the \BATSRUS\ code is called 8-wave scheme,
because the MHD equations can be rewritten in a form, where in addition
to the entropy, 2 slow, 2 Alfv\'en and 2 fast waves, there is an 8th 
wave associated with finite $\divB$. The momentum, 
energy , and induction 
equations 
 (\ref{eq:momentum_conserve})- (\ref{eq:energy_conserve}) 
and (\ref{eq:bfield_conserve}))
are modified to
\begin{eqnarray}
  \frac{\partial (\rho\bu)}{\partial t}+
      \nabla\cdot\left(\bu\rho\bu - \frac{\bB\bB}{\mu_0}\right) + 
      \nabla\left(p + \frac{\bB^2}{2\mu_0}\right)
      &=& -\frac{1}{\mu_0}(\divB)\bB                    \label{q-momentum8} \\
  \frac{\partial\bB}{\partial t} + \nabla\cdot(\bu\bB - \bB\bu)
      &=& -(\divB)\bu                    \label{q-induction8} \\
  \frac{\partial\varepsilon}{\partial t}+
      \nabla \cdot \left[ \bu \left( \frac12 \rho v^2 + 
      \frac{\gamma}{\gamma-1} p +
      \frac{B^2}{\mu_0} \right) - 
      \frac{\left(\bu\cdot\bB\right)\bB}{\mu_0}\right] 
      &=& -\frac{1}{\mu_0}(\divB)\bB\cdot\bu           \label{q-energy8}
\end{eqnarray}
The source terms on the right hand side are zero analytically, but
can be non-zero numerically due to truncation errors. The source
terms are added at every time step with some simple discretization.
The 8-wave scheme is very robust, and it seems to produce correct
results for the vast majority of MHD problems. A potential drawback,
however is, that it may not give correct jump conditions, since it
is not conservative.

\subsection{Diffusive Control \label{section:diffusive_control}}

Diffusive control (Linde and Malagoli, 2000) adds terms that diffuse
away the numerically generated $\divB$. The new source terms occur in 
the induction and energy conservation equations
\begin{eqnarray}
  \frac{\partial\bB}{\partial t} + \nabla \cdot(\bu\bB - \bB\bu)
      &=& D \nabla(\divB)                    \label{q-inductiondiff} \\
  \frac{\partial\varepsilon}{\partial t}+
      \nabla \cdot \left[ \bu \left( \frac12 \rho v^2 + 
      \frac{\gamma}{\gamma-1} p +
      \frac{B^2}{\mu_0} \right) - 
      \frac{\left(\bu\cdot\bB\right)\bB}{\mu_0}\right] 
      &=& \frac{1}{\mu_0} D \bB\cdot\nabla(\divB)             \label{q-energydiff}
\end{eqnarray}
Taking the divergence of the modified induction equation, we obtain
\begin{equation}
  \frac{\partial\divB}{\partial t}= \nabla\cdot D \nabla(\divB) 
\end{equation}
which is a simple diffusion equation for $\divB$ with a diffusion
coefficient $D$. The diffusion coefficient is set to a value that provides 
maximum diffusion of $\divB$ but still keeps the scheme numerically stable.
This requirement is met if
\begin{equation}
  D=\delta \frac{(\Delta x)^2}{\Delta t}
\end{equation}
with $\delta\le4/6$.  This value is marginally stable, so the value used
should be smaller than this.  We have found that numbers around 0.4 work
well.  More testing remains to be done to find the optimal parameter.

The diffusive source terms can be added to the 8-wave source terms, or used
by themselves (Sokolov). In the above form the diffusive source terms
are only conservative for the induction equation and for a uniform grid.
However we have some ideas on how to make the diffusive scheme fully 
conservative and these ideas will be implemented and tested soon.

\subsection{Projection Scheme \label{section:projection}}

The projection scheme eliminates the numerically generated $\divB$ in
every time step. The solution $\bB^*$ provided by the base scheme can
be written as the sum of a curl and a gradient
\begin{equation}
  \bB^*=\nabla\times\bA+\nabla\varphi         \label{q-curlplusgrad}
\end{equation}
where the physically meaningful and divergence free part is the first
$\nabla\times\bA$ term. By taking the divergence of the above equation,
we get a Poisson equation for the unknown scalar field $\varphi$:
\begin{equation}
  \divB^{*}=\nabla^2 \varphi                  \label{q-poisson}
\end{equation}
This equation can be solved with some iterative scheme. In \BATSRUS, the
BiCGSTAB scheme is used, which is a parallel Krylov-type method.
Once the solution is found (with some accuracy) the divergence free solution
can be easily obtained
\begin{equation}
  \bB^{n+1}=\bB^*-\nabla\varphi                \label{q-correction}
\end{equation}
It is interesting to note that the correction term above is very similar
to the diffusive source term in the modified induction equation 
(\ref{q-inductiondiff}) of the diffusive control scheme if we take 
$\varphi=D \divB$. This is a very crude solution, but  
for certain iterative schemes it would be the first iterate
of the Poisson equation.

\subsection{Constrained Transport Scheme \label{section:CT}}

The Constrained Transport (CT) scheme relies on a specific staggered 
discretization of the magnetic field and the electric field which maintains 
the divergence free property. {\bf It is important to note that 
THE CT SCHEME CAN ONLY BE USED IN THE TIME ACCURATE MODE}. 
It is also important to distinguish the fully conservative 
CT scheme (which follows Balsara and Spicer 1999 and was generalized to AMR by 
T\'oth and Roe, 2000) which is used in \BATSRUS, 
from the original non-conservative Hawley and Evans (1988) scheme.

A main feature of the CT scheme is the introduction of variables that
are stored on face centers and cell edges.  The previously described 
schemes all have  variables stored at the cell center.  For CT, the magnetic 
field is
stored at face centers (we use lower case $\bf b$ to distinguish from the
above cell centered $\bf B$) and we introduce the electric field $\bf E$
which is stored on cell edges.
The main steps of the algorithm are the following
\begin{itemize}\itemsep=0pt
\item Discretize $b^x$, $b^y$, $b^z$ at the $x$, $y$, $z$
             face centers of the cells, respectively
\item Calculate $E^x$, $E^y$, $E^z$ at cell edges by interpolation
      of fluxes provided by the base scheme
\item Message pass the electric field $\bE$ for consistency at refinement
      changes
\item Update $b^x$, $b^y$, $b^z$ with finite differencing $\nabla\times\bE$
\item Interpolate $\bb$ to the cell centered $\bB$
\end{itemize}

As long as $\divb^n=0$, this algorithm ensures that $\divb^{n+1}=0$ as well.
\BATSRUS\ uses the projection scheme to ensure that the initial condition
has exactly zero $\divb$. 
When ever a new run is begun, the code is restarted, or the $\divB$ control
method is changed to CT from some other scheme, it is recommended that 
a projection
be computed rather accurately to ensure that the initial condition for the CT
scheme is divergence free.  With this in mind, it is possible to switch to the 
CT scheme at any time during a run.

In a CT scheme the boundary conditions must be applied on the electric
field. Currently the code uses a simple $\bE=0$ condition at the inner
boundary. This should be improved in the future by using the electric
field provided by the ionosphere model.

Refinement and coarsening maintains $\divb=0$ using specially designed
restriction and prolongation operators. However, the blocks intersecting
the inner boundary should not be refined or coarsened while the CT scheme
is used, because the current algorithm cannot deal with the newly
created or removed cells while maintaining $\divb=0$. This might be
improved in the future.

\subsection{Choosing a Divergence B Control Scheme \label{section:choosing_divb_control}}

For steady state problems local time stepping is strongly 
advised. This excludes using the constrained transport scheme
for $\divB$ control. The projection scheme is quite expensive
for the AMR grid when the code is running on many processors.
The current implementation of the diffusive divergence control
seems to produce results which seem to be less accurate than the
results obtained with the other schemes. 
Therefore we recommend the 8-wave scheme as the most efficient method 
for obtaining a steady state solution.
With an improved implementation of the diffusive control it is
possible that the diffusive control in itself or combined with the
8-wave scheme will work better. 

If the steady state solution shows an excessive amount of $\divB$
(as compared with the current or the gradients of $\bB$) then 
adding diffusive control is likely to reduce $\divB$ by a factor
of about 10. The projection scheme can reduce the error even more, 
but it is costly.

For time accurate problems the constrained transport, the 8-wave
and the diffusive control need about the same CPU time, while the
projection scheme is still very expensive if \BATSRUS\ is run
with many blocks on many processors.
For most time accurate problems we find the 8-wave scheme to be the 
most robust discretization. If keeping $\divB=0$ has a high priority,
one can experiment with the constrained transport scheme. Note, however,
that the boundary conditions are not fully developed, and they are 
mostly appropriate for an Earth magnetosphere run without the ionosphere
model.

